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We present some advances, both from a theoretical and from a computa- 
tional point of view, on a quadratic vector equation (QVE) arising in Marko- 
vian Binary Trees. Concerning the theoretical advances, some irreducibility 
assumptions are relaxed, and the minimality of the solution of the QVE is 
expressed in terms of properties of the Jacobian of a suitable function. From 
the computational point of view, we elaborate on the Perron vector-based 
iteration proposed in p.]. In particular we provide a condition which ensures 
that the Perron iteration converges to the sought solution of the QVE. More- 
over we introduce a variant of the algorithm which consists in applying the 
Newton method instead of a fixed-point iteration. This method has the same 
convergence behaviour as the Perron iteration, since its convergence speed 
increases for close-to-critical problems. Moreover, unlike the Perron iteration, 
the method has a quadratic convergence. Finally, we show that it is possible 
to alter the bilinear form defining the QVE in several ways without changing 
the solution. This modification has an impact on convergence speed of the 
algorithms. 

Keywords: Markovian binary tree, branching process, Newton method. Perron vector, 
fixed-point iteration 

1 Introduction 

Markovian Binary Trees (MBTs) are a particular family of branching processes, which 
are used to model the growth of populations consisting of several types of individuals 
who evolve independently and may produce a variable number of offsprings during their 
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lifetime. MBTs have applications in biology, epidemiology and also in telecommunication 
systems. We refer the reader to [Si [3] for definitions, properties and applications. 

One important issue related to MBTs is the computation of the extinction probability 
of the population, which can be characterized as the minimal nonnegative solution 

with := {v eR^ : Vi > 0,i = 1, . . . , N}, 

of the quadratic vector equation (QVE for short) 

X = a + b{x, x), (1) 

where a G M^, b : M.^ x R^ M.^ is a vector- valued bilinear form, such that the vector 
e = (1, 1, ... , 1)"^ is always a solution of ([T]). For writing the entries of b in coordinates, 
we use the notation bijk ■= efb{ej,ek), where ee is the Ith vector of the canonical basis. 
With this choice, 

{b(x,y))i = ^bijkXjVk- 

In many papers the notation b{s,t) = B{s ® t), with B € M.^^^ and (8) denoting 
the Kronecker product, is used instead; one can see that the two representations are 
equivalent. We favor the former, since it highlights the symmetry features of the problem. 
We mention the fact that the functions obtained by fixing the first or the second argument 
of the bilinear form, i.e., b{y, •) and b{-,z) for suitable y,z & M+, are linear maps from 
to itself, and thus they can be represented by x matrices with nonnegative 
entries. 

The MET is called subcritical, supercritical or critical if the spectral radius p{R) of 
the matrix 

R:=b{e,-) + b{-,e) (2) 

is strictly less than one, strictly greater than one, or equal to one, respectively. 

Under the stated assumptions, one can prove the existence of a minimal nonnegative 
solution in the componentwise ordering. A proof using minimal hypotheses is presented 
in [4J. In the subcritical and critical cases the minimal nonnegative solution is the vector 
of all ones, while in the supercritical case x* < e, x* ^ e. Thus, only the supercritical 
case is of interest for the computation of x*. 

Moreover, in the following we shall focus on the case in which x* > 0. It is shown in 
[3] how to detect reliably the cases when this property does not hold, and reduce them 
by projection to problems of lower dimension with strictly positive minimal solution. 

Several iterative methods have been proposed and analyzed for computing the vector 
x*. In [2] the authors propose two fixed point iterations with linear convergence, called 
depth and order algorithms. Another linearly convergent algorithm, called thicknesses 
algorithm, is proposed in [3]. In [5] and in [6] two variants of Newton's method are 
proposed. A different algorithm, based on a Perron vector iteration, is proposed in [1]. 
This algorithm, unlike classical iterative methods, increases its convergence speed for 
close to critical problems. 
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In this paper we provide theoretical and computational advances concerning the QVE 
([T]). We first show that the matrix of ([2]) can be assumed irreducible, since if it 
were reducible, we may reduce the problem of solving ([TJ to the problem of solving 
QVEs of smaller dimension, whose associated matrix R is irreducible. Assuming that 
R is irreducible, we provide a new characterization of the minimal nonnegative solution 
X*, in terms of the properties of the Jacobian of the function F{x) = x — a — b{x,x), 
evaluated at x = x*. This property, which complements the results in [1], allows us 
to give a condition which ensures that the limit of the Perron vector-based iteration 
provides the sought solution x* of the quadratic vector equation. 

Moreover, we introduce a variant of the Perron vector-based iteration, which con- 
sists in applying the Newton method instead of a fixed-point iteration. This method is 
quadratically convergent, and has the same convergence behaviour as the Perron iter- 
ation, as its convergence speed increases for close-to-critical problems. The number of 
iterations needed by this variant is usually lower than the number of iterations needed 
by the original Perron iteration. However, due to the larger complexity of the single 
iteration step, the Newton-based method is generally slower than the Perron iteration, 
in terms of total computational time. 

Finally, we show that it is possible to alter the bilinear b{x, y) form definining the 
QVE in several ways without changing the solution. This modification has an impact on 
convergence speed: in most examples, making the wrong choice can double the number 
of iterations needed. We show that, at least on the experiments reported, the best results 
are given by a symmetrization of the original bilinear form. 

The paper is organized as follows. In Section [2] we recall classical algorithms based on 
fixed point iterations, while in Section [3] we recall the Perron-based iteration. In Section 
Hlwe discuss the case where the matrix i? of ([2]) is reducible, and we reduce the QVE to 
smaller size QVEs whose associated matrix R is irreducible. In Section [5] the minimality 
of the solution x* of the QVE is expressed in terms of properties of the Jacobian of the 
function F(x) at x = x*. This result is used in Section [U] to ensure that the limit of 
the Perron-based iteration provides the sought solution x* . The Newton version of the 
Perron-based iteration is proposed in Section [71 In Section [5] the choice of the bilinear 
form b{x,y) is discussed. The results of the numerical experiments are presented and 
discussed in Section O We draw conclusions in Section [TUJ 

2 Classical iterations 

Several iterative methods have been proposed and analyzed for computing the vector x*. 
In [2] the authors propose two iterations with linear convergence, called depth and order 
algorithms, which are defined respectively by the two linear equations 

(/ - b{-,Xk))xk+i = a, 
{I - b{xk, ■))xk+i = a. 

The thicknesses algorithm, still linearly convergent, is proposed in [3] and consists in 
alternating iterations of each of the two above methods. 
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In [5] the authors apply the Newton method to the map 



F{x) := X — a — h{x, x) 



(3) 



obtaining the iteration defined by 



(/ - h{xk, ■) - b{-,Xk))xk+i = a- b{xk, Xk) 



(4) 



which converges quadratically. Its convergence speed is usually much higher than that 
of the previous, linearly-convergent iterations. A modification of the Newton method, 
which increases slightly its convergence speed, has been proposed in [B]. 

All these methods have probabilistic interpretations, in that their k-th iterate Xk can 
be interpreted as the probability of extinction of the process restricted to a special subtree 
7fc. Each of them provides a sequence {xk}k of nonnegative vectors, with xq = (0, . . . , 0)"^, 
which converges monotonically to the minimal nonnegative solution x*. A common 
feature of all these methods is that their convergence speed slows down when the problem, 
while being supercritical, gets close to critical, i.e., the vector x* approaches the vector 
of all ones. This happens because the mean extinction time increases, and thus sampling 
larger and larger trees is needed to capture the behavior of the iteration. 

3 A Perron-vector-based iteration 

In [1], the authors propose an iterative scheme based on a different interpretation. Let 
us suppose for now that the nonnegative matrix R, as defined in ([2]), is irreducible — 
we discuss this assumption in Section SI Then, since irreducibility depends only on the 
zero pattern of the matrix, b{u, •) + b{-,v) is irreducible for each u,v a with strictly 
positive entries. 

If we set y = e — x, equation ([1]) becomes 



A special solution of ([5]) is y* = e — x* , where x* is the minimal nonnegative solution 
of ([1]). Notice that < y* < e. In the probability interpretation of Markovian Binary 
Trees, since x* represents the extinction probability, then y* = e — x* can be interpreted 
as survival probability. In particular, y* is the probability that a colony starting from a 
single individual in state i does not become extinct in a finite time. The three summands 
in the right-hand side of dSJ also admit an interesting probabilistic interpretation [T]. 
If we set Hy := b{-, e) + b{e — y, •), equation ([5]) becomes 



If Hy is nonnegative and irreducible (which happens for sure if y < e, in view of the 
irreducibility of R), then the Perron- Frobenius theorem implies that p{Hy*) = 1 and y* 
is the Perron vector of the matrix Hy*. 

This interpretation allows to design new algorithms for computing y* and x* . Applying 
a functional iteration directly to ([6]), or the Newton method, gives nothing new, since 




(5) 



y = Hyy. 



(6) 
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we just did a change of variable. However, if we define the map PV(M) as the Perron 
vector of a nonnegative irreducible matrix M, we may rewrite ([6]) as 

y = FY{Hy). (7) 

We may apply a fixed-point iteration to solve ((Tj), thus generating a sequence {yk}k of 
positive vectors such that the vector yk+i is the Perron vector of the matrix -ffy^., i.e., 

yk+i = PV{HyJ. (8) 

A suitable normalization of the Perron vector, consistent with the solution, is needed 
to obtain a well-posed iteration. An optimal normalization choice is suggested in [IJ. If 
we take w as the Perron vector of the nonnegative irreducible matrix , then we may 
normalize yk+i so that 

w^iVk+i - b{yk+i,e) - 6(e, y^+i) + b{yk+i,yk+i)) = 0, (9) 

i.e., we impose that the residual of Q for y = yt+i is orthogonal to w. With this 
choice, one can prove [T] that the convergence speed of the sequence {yk}k defined in 
([8]), with the normalization condition ([9]), is linear with a small convergence factor for 
close-to-critical problems, and tends to super linear as the considered problems approach 
criticality. Thus, although the convergence of this method is linear, surprisingly its speed 
increases as the problem gets close to critical, unlike the classical iterations. 



4 Dealing with reducible R 

The following result shows that when R is reducible we can reduce a QVE to two smaller- 
dimension problems to be solved successively with a kind of back-substitution. Therefore, 
for the solution of a generic QVE, we only need to apply the Perron iteration to the case 
in which R is irreducible. 



Theorem 1. Suppose that, for a QVE ([T]) with x* > 0, we have 

R 



Rii Ri2 
R22 



where Rn is M x M and R22 is {N - M) x {N - M). Let 



Xi 
X2 



be partitioned accordingly. Let 



P := [I 



M 







MxiN-M)\ 1 



Q:=[0 



{N-M)xM 



L 



N-M \ , 



be the orthogonal projection on the first M and last N — M components respectively. Let 
us define the bilinear form on 



b2iu,v) ■.= Qb{Q^u,Q^v). 
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Moreover, for each y G M^^^^, let us define 

Ty ■.=Im - Pb{;Q^y)P^ - Ph{Q^y, ■)P^, ay ■.=Ty-\ai + Pb{Q^ y , y)) , 
and the bilinear form on M^^ 

by{u,v) := Ty-^Pb{P^u,P'^v). 

Then, 

1. X solves ^ if and only ifTx2 is nonsingular, and its block components X2 and xi 
solve respectively the two quadratic vector equations 

X2 = a2 + b2{x2,X2) (10) 

and 

xi = a^2 + 6x.2(xi,xi). (11) 

2. If X* is the minimal solution to then and are the minimal solution to 
(fTTj) and ([To]) respectively. 



Proof. First notice that since P and Q are projections on complementary subspaces, ([T]) 
holds if and only if it holds when projected on both, i.e., 

xi =ai + Pb{x,x), (12a) 
X2 =02 + Qb{x,x). (12b) 

Since Rij = Ylk=ii^ijk ~^bikj) , it follows from the block structure of R (and from bijk > 0) 
that bijk = whenever i > M and either j < M or k < M. This implies that the second 
block row of b{u, v) depends only on the second block rows of u and v. We can write 
this formally as Qb{u,v) = Qb{Q^Qu,Q^Qv). Then, ()12bp is equivalent to ^TO\i . 

By exploiting bilinearity and the fact that P^P + Q^Q = In, we can rewrite ()12aP as 

XI = ai + P6(P^xi, P^xi) + P6(P^xi, Q^xa) + Pb{Q^X2, P^xi) + Pb{Q^X2, Q^xs), 



or 

T^^xi = oi + Pb{Q^X2, 0^X2) + Pb{P^xi,P'^xi). 

Since the right-hand side is nonnegative and xi is positive, the Z-matrix T^^ is an M- 
matrix. Therefore, we may invert it to get ()lip . The steps in the above proof can be 
reversed provided is nonsingular, thus the reverse implication holds as well. 

Let us now prove the second part of the theorem. Equation (jlOp admits a minimal 
solution due to the general existence theorem (since it admits at least a solution); suppose 
it is X2 7^ X2; then, by minimality, X2 < X2. The matrix > T^* is an M- matrix, thus 
T~} > . Therefore, < o^;* and bx2 < b^*. We have 

X\ = a^x* + bx* {xi,Xi) > 0x2 + (^1)^1)) 
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thus the equation (jlip has a supersolution, and this imphes that it has a solution by [U 
Lemma 5]. Let xi be its minimal solution; then, x is a solution to ([1]) by the first part 
of this theorem, but this is in contradiction with the minimality of x*, since X2 ^ X2- 
Therefore X2 is the minimal solution to (jlOp . If ()lip admitted a solution xi ^ xj, then 
by the first part of the theorem 

would be a solution to ([1]) , and this again contradicts the minimality of x* . □ 

Let := I — b{x,-) — b{-,x) be the Jacobian of the map F{x) defined in ([3|) (see 
[5111]). Notice that if x > 0, then F^ has the same positivity pattern as R, and thus is 
irreducible whenever R is. Moreover, when x = e is a solution to ([1]), then the all-ones 
vectors of suitable dimension solve (|10p and (jlip , thus the Perron vector-based iteration 
can be applied to the reduced problems as well. 

5 An alternative characterization of minimality 

The following theorem provides a practical criterion to check the minimality of a solution. 

Theorem 2. Let x > be a solution of ([1]) and assume that R is irreducible. Then, F^ 
is an M-matrix if and only if x is minimal. 

Proof. The implication (x* minimal) =^ (-F^., is an M-matrix) has been proved in [4J. 
We prove the reverse here. The proof is split in two different arguments, according to 
whether F^ is a singular or nonsingular M-matrix. 

Let F^ be a nonsingular M-matrix, and let x be another nonnegative solution; we need 
to prove that x — x > 0. Prom the Taylor expansion of F{x) (and the fact that F" < 0) 
we have 

1 

= F(x) = F(x) + F^{x -x) + ^Fx{x - x, x - x) < F^{x - x), 

that is, F^{x — x) > 0. It suffices to multiply by {F!^)~^ > to get x — x > 0, as needed. 

Let now F^ be a singular M-matrix. Suppose that x is not minimal, and x* ^ x is the 
minimal solution to ([1]). Then, F^, ^ is a (singular or nonsingular) M-matrix, by the 
reverse implication of this theorem. Thus, by the properties of M-matrices, F^ must be 
a nonsingular M-matrix, which is a contradiction. □ 

Notice that this characterization of minimality allows to deduce easily the fact, claimed 
above, that the solution e is minimal only in the subcritical and critical cases. 

6 On the limit of the Perron iteration 

The following result shows that, under reasonable assumptions, the limit of the Perron 
vector-based iteration is the minimal solution of ([T]). 
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Theorem 3. Suppose that R is irreducible, and that x* > 0. Suppose that the Perron 
iteration ([8]), with normalizing condition ([9]), converges to a vector y* such that y* < e. 
Then, x = e — y* is the minimal solution of ([1]). 

Proof. Let us first prove that the spectral radius of Hy* is 1 . The iterates of the Perron 
iteration satisfy 

Xk+iVk+i = Hyi^yk+i, (13a) 
w'^iVk+i - = 0. (13b) 

By passing (fT3|) to the limit, we get 

X*y* = Hy.y\ (14a) 
uF{y* - Hy,y*) = 0. (14b) 

Notice that A* is well-defined, as it may be defined as the common ratio between the 
components of Hy*y* and those of y* . We left-multiply (|14ap by uF to get uF {\*y* — 
Hy*y*) = 0, which, compared to ()14bp . tells us that A* = 1. In particular, this implies 
that X = e — y* is a solution of ([T|), as we may verify directly by back-substitution. 

Moreover, p{Hy*) = 1, and thus / — Hy* = I — b{e — y* ,■) — b{-, e) is a singular M- 
matrix. Thus the Z-matrix F!^ = I — b{e — y*, •) — 6(-, e — y*) > I — b{e — y*, •) — 5(-, e) 
is an M-matrix, too. By Theorem [21 this implies that x = e — y* is minimal. 

□ 



7 The Perron-Newton method 

We may also apply Newton's method for the solution of ([7]). 

We first recall the following result from [1], which provides an explicit form for the 
Jacobian of the iteration map G{y) defining the iteration ([8]) with the normalization ([9]), 
i.e., 

G{y) := the Perron vector of Hy, normalized s.t. w'^ {G{y) — HQ(^y-jG{y)) = 0. 

Theorem 4. Let y be such that Hy is nonnegative and irreducible. Let u = G{y), and 
let V be such that v'^Hy = Xv^ , where A = p{Hy). Then the Jacobian of the map G at y 

T T 
JGy =(^L- ^) {Hy - A/)t (^L - b{;u), (15) 

where 

:= w'^{J — b{e — u,-) — b{-, e — u)) 
and the symbol ^ denotes the Moore-Penrose pseudo-inverse. 
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input : the bilinear form b (note that a is not necessary — in fact it can be 

deduced from e = a + 6(e, e)) 
input : the normahzation vector w > {a good choice is taking the Perron vector 

of i?^, see pP) 

y ^ e; 

while a suitable stopping criterion is not satisfied do 
u ^ G{y); 

J ^ JGy (computed using p^): 

y ^y-{i - jy^iv - u); 

end 

if < y < e then 
I X ^ e-y; 
else 

I (error: no convergence); 
end 

Algorithm 1: The Perron-Newton algorithm 



With the aid of this formula, we may define the Perron-Newton method for the solu- 
tion of (dl) as in Algorithm [TJ 

A step of Newton's method basically requires a step of the Perron vector-based fixed- 
point iteration associated with ([7]), followed by the computation of a Moore-Penrose 
pseudoinverse and the solution of a linear system. Thus its cost is larger than, but still 
comparable to, the cost of a step of the Perron vector-based functional iteration. This is 
compensated by the fact that the Newton method has quadratic convergence, and thus 
requires less iterations. 

The convergence properties of the Perron-Newton method for close-to-critical prob- 
lems are similar to those of the Perron vector-based functional iteration. For close- 
to-critical problems one has x* ~ e, therefore p{JGy*) ~ 0. Hence, by the Newton- 
Kantorovich theorem [7] there is convergence for sufficiently close-to-critical problems. 
The proof of Theorem [3] can be easily adapted to show that the limit point must corre- 
spond to the minimal solution of ([1]) if < < e. Moreover, since p{JGy*) ~ 0, the 
matrix to invert is well-conditioned and y — G{y) has a simple zero. 

8 On the choice of the bilinear form h 

Equation ([1]), and thus its solution, depend only on the quadratic form 6(t, t) := B{t®t); 
however, there are different ways to extend it to a (nonnecessarily symmetric) bilinear 
form b{s,t). Namely, for each i and each j ^ k, we may alter simultaneously bijk 
and bikj, as long as their sum remains invariant, and they both remain positive. For 
example, we may switch the two terms in every such pair, obtaining the bilinear form 
b'^{s,t) := b{t,s). 

Some of the solution algorithms depend essentially on the choice of the bilinear ex- 
tension: for instance, the depth and order algorithms. It is easy to see that the depth 
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algorithm applied to b'^ coincides with order applied to 6, and vice versa. Instead, in 
the classical Newton's method the bilinear form appears only in the expressions 
bixk, ■) + b{-,Xk) and b(xk,Xk), which are unaffected by this change. Thus the classical 
Newton method stays the same no matter which bilinear extension we choose. 

On the other hand, one can see that the Perron-vector based functional iteration and 
its Newton based version do depend on the bilinear extension, and their convergence 
speed is affected by this choice. The expression of the bilinear form ultimately reflects a 
modeling aspect of the problem. While in the original definition of a branching process 
an individual splits into two new ones in two different states, it is often convenient to 
identify one as the "mother" and one as the "child" , even if this distinction is artificial. 
In fact, we can safely redefine who is the mother and who is the child, as long as we do 
not change the total probability that an individual in state i originates two offsprings in 
states j and k. This corresponds exactly to changing the bilinear form b in the described 
way. 

Among the possibilities for the modifications of b, we list the following. 
Transposition b'^{s,t) ■.= b{t,s) 
Symmetrization b^{s,t) := \ {b{s,t) + 6^(s,t)) 
Desymmetrization 1 

bijk + bikj if j < k 
bijk if j = k 

ifj>k 

Desymmetrization 2 

bijk + bikj ifj>k 
bijk ifj = k 

ifj<k 

In the following section, we report numerical experiments performed with the above 
bilinear extensions and compare the computational times. We do not have a definitive 
result on which choice gives the best convergence: as is the case with the depth and 
order algorithms, the optimal bilinear extension may vary in different instances of the 
problem. 

9 Numerical experiments 

We performed numerical experiments to assess the speed of the proposed methods. The 
tests were performed on a laptop (Intel Pentium M 735 1.70Ghz) with Matlab R2010a 
and considered two sample parameter-dependent problems. 

PI a small-size Markovian binary tree with branches of varying length, described in [5l 
Example 1]. It is an MBT of size = 9 depending on a parameter A, which is 
critical for A ~ 0.85 and supercritical for larger values of A. 



{b^'^k ■■= (&^) 1 = 
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Figure 1: CPU time vs. parameter A for PI — lower=better 



P2 a random-generated MBT of larger size {N = 100). It is created by generating a 
random bilinear form b, choosing a suitable a so that a + b{e, e) = Ke for some 
K, and then scaling both a and b in order to eliminate K. We report the Matlab 
code used for its generation in Algorithm [2j Larger choices of the parameter A 

input : the size of the MBT and a parameter A > 

e=ones(N,l); 

rand ('state ',0) ; 

b=rand(N,N*N); 

K=max(b*kron(e,e))-|-$\lambda$; 

a=K*e— b*kron(e,e) ; 

a=a/K; 

b=b/K; 

Algorithm 2: Generating a random MBT 

increase the values of a, i.e., the probability of immediate death, and thus enlarge 
the extinction probability making the process closer to critical. With A^ = 100, 
the process is critical for A ~ 4920. 

Figure [T] shows a plot of the computational times for classical Newton and the two 
Perron vector-based methods for different values of the parameter A. Depth, order and 
thicknesses are not reported in the graph as they are much slower than these methods, 
as also shown by the experiments in f5]. While in close-to-critical cases the time for CN 
has a spike, the ones for PN and PI seems to decrease. While having in theory worse 
convergence properties, the Perron iteration is faster than the Perron Newton method: 
the additional overhead of the pseudoinverse and of the computation of both left and 
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Figure 2: CPU time vs. parameter A for P2 
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Figure 3: Number of steps needed for the Perron iteration for PI with several variants 
of the bilinear form 
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Figm'e 4: Number of steps needed for the Perron iteration for P2/ with several variants 
of the bilinear form 



right dominant eigenvector do not make up for the increased convergence rate. 

Figure [2] shows the corresponding plot for the larger problem P2. We point out that 
two different methods were used to compute the Perron vectors in the two problems. 
For P2, we use eigs, which is based on an Arnoldi method [8j. On the other hand, for 
PI, due to the really small size of the problem, it is faster to compute a full eigenvector 
basis with eig and then select the Perron vector. 

With this choice, both the Perron iteration and Perron-Newton method are faster 
than the classical Newton method on this larger-size problem, in the close-to-critical 
region. 

Figure [3] reports the number of iteration (which essentially grows as the CPU time) 
for the Perron iteration on PI with several alternative bilinear forms equivalent to b. We 
see that among the two possible "branch switches" , in this example the iteration with b 
converges faster than the one with b'^. Clearly this cannot be a general result: due to 
the involutory nature of this transposition operation, if we started with b := 6^, then the 
faster choice would have been b^ = (6"^)^ = b. Thus we cannot infer a rule for telling 
which of the two is preferable. Similarly, it is impossible to do a proper comparison 
among b^^ and B^^. On the other hand, an interesting result is that the performance 
of the iteration with b^ seems to be on par with the better of the two. 

Figure H] reports the same comparison for the problem P2. The results are less pro- 
nounced than on the previous example: since the entries of the bilinear form b are gen- 
erated randomly, the difference between the "left" and "right" branches of the binary 
tree should be less marked than in PI, where the two directions are willingly unbalanced. 
Nevertheless, the symmetrized bilinear form consistently yields slightly lower iteration 
counts. 
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Therefore, based on these results, we suggest to apply the Perron iteration and Newton 
methods on the symmetrized bilinear form instead of the original one. 

10 Conclusions 

In this paper we presented several possible implementation variants of the Perron vector- 
based iteration introduced in [I]. A Newton method based on the same formulation of 
the problem is slightly less effective than the original iteration, although it maintains the 
same good convergence properties for close-to-critical problem. Moreover, we highlight 
the fact that there is a family of possible modifications to the bilinear form b that alter 
the form of solution algorithms, but not the original equation ([1]) and its solution. One 
of these modifications, the symmetrization, seems to achieve better results than the 
original formulation of the numerical algorithms. 

Moreover, we present a couple of theoretical results on quadratic vector equations that 
show how to ensure that the obtained solution is the desired one, and how to deal with 
the problems in which an irreducibility assumption is not satisfied. 
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